clear; clc; close all;

% 参数设置（同Fig5a）
q = 0.75;     
a = 2;       
P = 2;
Q = 1;
x0 = 2;
y0 = -2;

b_min = -5;
b_max = 5;
M_b = 100;               % MLE计算可采样点较少，避免过慢
b_list = linspace(b_min, b_max, M_b);

N_total = 1500;          % 总迭代步数
N_transient = 500;       % 丢弃暂态步数

MLE_list = zeros(1, M_b);

for idx = 1:M_b
    b = b_list(idx);
    [x, y] = SCLMM(q, a, b, P, Q, x0, y0, N_total);
    
    % 丢弃暂态后计算MLE
    X = [x(N_transient+1:end); y(N_transient+1:end)];
    
    % 计算最大李雅普诺夫指数（Wolf算法）
    MLE_list(idx) = wolfMLE(X, 10);
    
    fprintf('b=%.3f, MLE=%.5f\n', b, MLE_list(idx));
end

% 绘图
figure('Position',[200 200 900 600]);
plot(b_list, MLE_list, 'r-', 'LineWidth', 1.5);
xlabel('b', 'FontSize', 14);
ylabel('Maximum Lyapunov Exponent (MLE)', 'FontSize', 14);
title(sprintf('Fig.5(b) 2D-SCLMM 最大李雅普诺夫指数 (q=%.2f)', q), 'FontSize', 16);
grid on;
box on;

%% 分数阶SCLMM函数
function [x, y] = SCLMM(q, a, b, P, Q, x0, y0, N)
    x = zeros(1, N);
    y = zeros(1, N);
    x(1) = x0;
    y(1) = y0;

    w = zeros(1, N);
    for k = 1:N
        w(k) = exp(gammaln(k + q - 1) - gammaln(k));
    end

    delta_x = zeros(1, N);
    delta_y = zeros(1, N);

    for n = 2:N
        delta_x(n-1) = a * sin(P * y(n-1)) * sin(P * x(n-1)) - x(n-1);
        delta_y(n-1) = b * cos(Q * x(n-1)) * cos(Q * y(n-1)) - y(n-1);

        sum_x = 0;
        sum_y = 0;
        for i = 1:n-1
            sum_x = sum_x + w(n - i) * delta_x(i);
            sum_y = sum_y + w(n - i) * delta_y(i);
        end

        x(n) = x0 + sum_x / gamma(q);
        y(n) = y0 + sum_y / gamma(q);
    end
end

